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ABSTRACT 

The amplitude and phase modulation observed in a significant fraction of the 
RR Lyrae variables - the Blazhko effect - represents a long-standing enigma in stellar 
pulsation theory. No satisfactory e xplanation for t he Blazhko effect has been proposed 
so far. In this paper we focus on the lStotherd (|2006l ) idea, in which modulation is caused 
by changes in the structure of the outer convective zone, caused by a quasi-periodically 
changing magnetic field. However, up to this date no quantitative estimates were 
made to investigate whether such a mechanism can be operational and whether it is 
capable of reproducing the light variation we observe in Blazhko variables. We address 
the latter problem. We use a simplified model, in which the variation of turbulent 
convection is introduced into the non-linear hydrodynamic models in an ad hoc way, 
neglecting interaction with the magnetic field. We study the light curve variation 
through the modulation cycle and properties of the resulting frequency spectra. Our 
results are compared with Kepler observations of RR Lyr. We find that reproducing 
the light curve variation, as is observed in RR Lyr, requires a huge modulation of 
the mixing length, of the order of ±50 per cent, on a relatively short time-scale of 
less than 40 days. Even then, we are not able to reproduce neither all the observed 
relations between modulation components present in the frequency spectrum, nor the 
relations between Fourier parameters describing the shape of the instantaneous light 
curves. 
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1 INTRODUCTION 

The periodic amplitude and phase modulation observed 
in many RR Lyrae variables - the Blazhko effect - is 
one of the most important, still unsolved problems in 
classical pulsation theory. In the last few years, thanks 
to extensive ground-based obse r vation campaigns (e.g., 
iKolenberg fc Guggenbergerl 120071: lJurcsik et al.l 120091 ) and 



satellite missions, CoRo T jchadid et al.ll2010h and Keple 
(e.g.. iBenko et al]|2010l ). we finally obtained nearly contin- 
uous data, allowing for detailed study of the light varia- 
tion associated with the Blazhko cycle. One of the most 
exciting new findings i s the period doubling phenomenon 
<|Kolenberg et alj|2010ah . never detected from ground-based 
data. It occurs at some phases of the Blazhko cycle and 
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manifests in an alternating shape of the light variation in 
consecutive pulsation cycles (see also lSzabo et al.ll2010l ). It 
was not detected in any non-modulated RR Lyrae star so 
far. Satellite data allow to study the light curve changes over 
the Blazhko cycle in detail. Cl early, consecutive Blazh ko cy- 
cles are not exactly repetitive (|Kolenberg et al1l201lT ). Also, 
the period of the Blazhko modulation differs from cycle to 
cycle. In the frequency spectra, the Blazhko phenomenon 
manifests itself as equidistantly spaced multiplets around 
main pulsation frequency and its harmonics. Triplets and 
quint uplets are detected in ground-based obse rvations (see 
e.g.. lJurcsik et al.|[20oi : IKolenberg et all 120091 ). In satellite 
data, not only triplets and quintuplets are visible, but also 
tenth-order side peaks (|Chadid et al.ll2010l ). 

All the newly discovered features of the Blazhko ef- 
fect put stringent constraints on the models proposed 
to explain this longstanding enigma. In fact, the two 
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most elaborated models, th e Magnetic Obli que Rota- 
tor/Pulsator model (MORP, IShibahashil I2000T ). and the 
Non-radial Resonant Rotator/Pulsator model (NRRP, e.g. , 
Van Hoolst et al.1 ll99Sl; iNowakowski fc Dziembowskil l200ll: 



Dziembowski fc Mizerskil 12004 ) are ruled out by observa- 
tions for several reasons. First, strong dipole magnetic fields, 
necessary for the MORP model to wo r k, were not detected in 
RR L yrae stars l|Chadid et al.l |2004| ; iKolenberg fc Bagnulol 
l2009h . Secondly, both the MORP and NRRP models pre- 
dict that the observed light curve modulation should mani- 
fest through specific features in the frequency spectra. The 
dipole geometry of the magnetic field, considered in the 
MORP model, leads to a quintuplet structure in the fre- 
quency spectra. A more complicated field geometry, that 
could possibly lead to higher order modulation components 
was not considered so far. Within the NRRP model it is 
shown that the modes of / = 1 are most easily excited, giv- 
ing rise to a triplet structure in the frequency spectra. The 
excitation of higher order nonradial modes is less likely. In 
addition, due to cancellation effects, the expected ampli- 
tudes would be very low for high-order modes. To the con- 
trary, in satellite dat a we clearly d et ect h igher-order mod- 
ulation components. IChadid et al.l l|2010t) detected mod- 
ulation components up to the tenth order. In both the 
MORP and NRRP models, strong asymmetries in the am- 
plitudes of the side peaks, as observed in the majority 
of Blazhko stars, are not trivial to reproduce. Consider- 
ing the triplets, in about 75 per cent of the Blazhko vari- 
ables, the higher frequency side peak has a higher am- 
plitu de than the lower frequency side peak ijAlcock et al.l 
2003). Third, both models propose mechanisms that imply 
a clockwork regularity, predicting robust modulation peri- 
ods, e.g., equal to the rotation period of the star, and re- 
peatable Blazhko cycles. Neither is seen in the data. The 
Blazhko vari ation can change consi derably even from cy- 
cle to cycle ()Kolenberg et all l201ll ) . In several stars two 
modulation per i ods a r e apparent (e.g., Detre fc Szeidl|[l973l ; 
lLaCluvze et all 120041 : ISzczygiel fc Fabrvckvl 120071 ). Hence, 
the connection between the period of the Blazhko modula- 
tion and the rotation period is highly questionable. Also, in 
several cases consecutive Blazhko cycles differ significantly, 
which manifests, e.g., in systematic growth or decay of the 
modulation amplitude. 

We also mention the recent s t udies of the Blazhko phe- 
nomen on by iBuchler fc Kollathl l|201ll ) and ljurcsik et alJ 
d201ll). Using the amp lit ude equation form a lism (see e.g., 
IBuchler fc Goupi]||l984l ). IBuchler fc Kollathl (|201ll ) showed 
that the high ord er half-integer resonance, proposed by 
ISzabo et al] (|2010h to explain the period doubling, can also 
cause the phase and amplitude modulation as observed in 
Blazhko variables. The result is very exciting and a more de- 
tailed analysis and confirmation t hrough the hydrody namic 
modelling would be of great value, [jurcsik et al.l l|201ll ) stud- 
ied the Blazhko variables in the globular cluster M5. They 
showed that the Blazhko stars tend to be situated on the 
zero-age horizontal branch and at the blue edge of the fun- 
damental mode instability strip. They speculate that this 
location hints that the Blazhko effect may be connected to 
the mode switch from the fundamental mode to the first 
overtone pulsation during the evolution. 

V ery recently, an idea proposed by IStothersI l| 200(1 
l2010h . which connects the Blazhko effect with variable mag- 



netic stellar cycles, has gained popularity. In this idea, the 
Blazhko modulation is connected to the cyclic strengthen- 
ing and weakening of turbulent convection in the outer stel- 
lar layers, caused by the postulated transient magnetic field. 
The field decays cyclically and is subsequently built up anew 
by the turbulent/rotational dynamo. Details of the process, 
particularly the interaction between magnetic filed, turbu- 
lent convection and pulsation were not discussed. Also, nu- 
merical estimates, e.g., on the expected strength of modula- 
tion or the necessary strength of the variable magnetic field, 
were not presented. For a critical analy sis of the Stothers 
idea we refer the reader to iKovac 1 (12009D. 

In this paper we use our pulsation hydrocodes to in- 
vestigate whether variable turbulent convection may cause 
such a light variation as is associated with the Blazhko effect. 
The variation of turbulent convection is put into the hydro- 
dynamical models in an ad hoc way - interaction with the 
magnetic field is neglected. Currently, there are no models 
capable of reproducing a variable magnetic field due to dy- 
namo mechanisms and of describing its dynamical coupling 
with turbulent convection and highly non-linear pulsation, 
features we deal with in the case of RR Lyrae variables. As a 
consequence, our model is strongly simplified and our results 
will serve for qualitative comparison with the observations. 
Having this in mind, we state the premise of our paper. If we 
can reproduce the most important observational constraints, 
the Stothers idea is certainly worth further, more detailed 
investigation. But if we do not succeed it needs revision. 

The structure of the paper is the following. In Sec- 
tion [5] we present the details of our numerical model testing 
the idea proposed by Stothers (2006) . In Section [3] we de- 
scribe the properties of our models and compare our results 
with observations, focusing on overall properties of Blazhko 
variables and on recent Kepler observatio ns of the famous 
Blazh ko variable, RR Lyr (KIC 7198959) (|Kolenberg et al.1 
2011). Some additional models are discussed in Section 2] 
and conclusions are presented in Section [5] 



NUMERICAL METHODS 

The basic tool in our modelling is the Warsaw non-linear 



convective pulsation hydrocode i|Smolec fc M oskalik 2 008a! ) ■ 
which we briefly describe in Section \2. II The code is slightly 
modified in order to model the effects of variable turbulent 
convection. Details, as well as the modelling procedure, are 
outlined in Section 12.21 In Section 12.31 we present the com- 
puted sequences of models to be analysed in this paper. 

2.1 Non-linear convective hydrocodes 

In all our computations we u se pulsation codes described 
by ISmolec fc Moskalik! i|2008al ). These are a static model- 
builder, linear non-adiabatic code and a direct time integra- 
tion non-li near hyd r ocode . For convective energy transfer 
we use the iKuhfufj { 19861 ) convection model reformulated 
for the use in stellar pulsation codes. Radiation is described 
in the diffusion approximation. The codes use a simple La- 
grangian mesh. 

For an extensive description and details of numerical 
implem entation we refer the reader to ISmolec fc Moskalik! 
|2008al ). Here we note that the model equations contain 
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eight order-of-unity scaling parameters that describe the 
turbulent convection model. These are mixing- length, a, and 
parameters multiplying the turbulent fluxes and terms that 
drive/damp the turbulent energy, a p (turbulent pressure), 
a m (eddy-viscous dissipation), a c (convective heat flux), a t 
(kinetic turbulent energy flux), a s (buoyant driving), aa 
(turbulent dissipation) and 7 r (radiative cooling). Theory 
provides no guidance for their values. However, some stan- 
dard values are in use, which result from comparison of 
a static, time-independent version of the model with the 
standard mixing-length theory ( see IWuchterl fc Feuchtingerl 
ll998l ; ISrnolec k, Mos kalik 2008a]). In practice, values of these 
parameters should be such that models satisfy as many ob- 
servational constraints as possible. 

We also stress that we use the original iKuhfufj (1 19861 ) 
prescription, in which essential buoyancy effects are included 
in convectively stable regions of the model (negative buoy- 
ancy). We note that the neglect of negative buoyancy (as 
is do ne, e.g., in the Florida-Budapest code, iKollath et al.l 
2002) can lead to unphysical effects in the computed mod- 
els. For an ex tensive discussion on the su bj ect we refer th e 
reader to, e.g JSmolec fc Moskalik! (|2008rJ ) or lSmoled (2009). 
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Figure 1. Possible shapes of mixing- length modulation consid- 
ered in this study. 



2.2 Exploring Stothers' idea 

IStothersl (2006) proposed that the Blazhko modulation is 
connected to the postulated cyclic strengthening and weak- 
ening of turbulent convection in the outer stellar layers, 
caused by a transient magnetic field. Turbulent convection 
becomes more vigorous during the decay of the magnetic 
field and it is quenched as the magnetic field builds up anew. 
Details of the underlying processes, particularly the interac- 
tion between magnetic field, turbulent convection and pul- 
sation were not elaborated by Stothers. 

Our code neglects the effects of magnetic fields and as 
such is not suitable for modelling the dynamical coupling 
of pulsation and turbulent convection on th e one hand, and 
the transient magnetic fields postulated bv lStothersI (|2006T l 
on the other hand. We only assume that the strength of the 
turbulent convection varies in time and do not elaborate on 
the physical mechanism behind such changes. The strength 
of turbulent convection is changed by cyclic variation of one 
of the Q-parameters entering our model (see Section l2.ip . 
The mixing length parameter, a, is our first choice, as it 
regulates the overall strength of convection and affects all 
the turbulent quantities enterin g the model (see equations in 
e.g. , ISmolec fc Moskalikll2008ah . We assume that the mixing 
length is either a sinusoidal or a triangular function of time, 
as illustrated in Fig. [1] 

The initial steps in our modelling correspond to a stan- 
dard hydrodynamic model computation, without modula- 
tion of the turbulent convection (constant mixing- length). 
First, we construct a static equilibrium model. Next, we 
compute the non-linear model. The static model is per- 
turbed with a scaled linear velocity eigenfunction followed 
by time evolution of the model. The model integration is 
stopped, when the full-amplitude single-periodic pulsation 
(the limit cycle) is reached. Then we start to modulate the 
turbulent convection. The non-linear model integration is 
restarted with the mixing-length parameter varying in time 
according to the chosen functional form (see Fig. [T]). The 
model is integrated for several Blazhko cycles. After a few 



initial cycles, the consecutive ones are almost indistinguish- 
able from one another, indicating that the Blazhko limit 
cycle is reached. Then, the model integration is stopped and 
the resulting light variation is subjected to a detailed anal- 
ysis (Section EI. 

In our computations we start to modulate the turbu- 
lent convection at the phase of maximum radius - the initial 
phase of modulation relative to the unperturbed mono-mode 
pulsation is equal to zero. To check whether the choice of the 
initial phase affects the results, we have computed several 
additional models with different initial phases: 0.25, 0.50 
and 0.75, and starting the modulation of turbulent convec- 
tion with an either increasing or decreasing mixing-length 
parameter. The final pulsation state (the Blazhko limit cy- 
cle) is insensitive to the initial phase. The computed light 
curves are almost indistinguishable for the models with dif- 
ferent initial phases (the trajectories overlap in the Fourier 
parameter plots, Figs. EHSJ). The amplitudes of the peaks in 
the frequency spectra (see Section f3.4p are also identical to 
within a fraction of a per cent. 

In Fig. [2] we present the light variation for a typical 
model over slightly more than one Blazhko cycle. The mod- 
ulation of the pulsation amplitude is clearly visible. In the 
lower panel we present the close-up of maximum amplitude 
phases at which the period doubling occurs0 



1 We note that our model predicts strictly periodic Blazhko cy- 
cles, which result from our assumption of a strictly periodic modu- 
lation of the mixing-length. This is not a necessary feature of the 
model. We could easily modulate the mixing-length in a quasi- 
periodic manner (as proposed in Stothers' paper), obtaining a 
quasi-periodic modulation. We have chosen periodic modulation 
for simplicity. 

2 See also the animated gif available in the electronic version 
of the Journal. It shows the light curve variation through the 
Blazhko cycle, including the period doubling phenomenon, as well 
as the variation of the Fourier parameters. 
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Figure 2. Bolometric magnitude versus time for a particular 
model of set M6 (T^ = 6800 K) . In the lower panel we show 
the close-up of phases at which period doubling occurs. 



Table 1. Physical parameters of the computed sequences of initial 
models. The chemical composition of all our models is the same, 
X = 0.76 and Z = 0.001. Effective temperatures are in the range 
of 6300 - 7000 K. 



Set 


M [Mq] 


L[L ] 


A 


0.65 


50.0 


13 


0.65 


40.0 


C 


0.65 


60.0 


D 


0.65 


70.0 



2.3 Computed model sequences 

In our modelling we focus on the fundamental mode mod- 
els only, as most Blazhko variables are fundamental mode 
pulsators. The chemical composition of all our models is 
the same, X = 0.76 and Z — 0.001, which is typical for 
RR Lyrae stars. We note that the metallicities of the Blazhko 
variables do not dif fer from thos e of the non-modulated 
RR Lyrae pulsators |SmolecJ 20051). In al l model computa- 
tions we use OP opacities llSeatonl I2005T) at low tempera - 
tures supplemented with the lAlexander fc Ferguson! (|l994j ) 
opacit y data. Opaciti e s wer e computed for the solar mix- 
ture of lAsplund et all (|2004t ). Other physical parameters of 
the initial, non-modulated models (masses and luminosities) 
are collected in Table [1] Parameters of set A, which is the 
most explor ed set in this study, ar e close to those quoted 
for RR Lyr (jKolenberg et alj|2010bh . In sets B, C and D the 
luminosity is varied, which allows to study the effects of a 
different M/L-ratio. 

The convective parameters of the initial non-modulated 
models are collected in Table [2] We adopt only one set 
of convectiv e parameters, very simila r to t he one w e hav e 
adopted in iBaranowski et all (|2009h and ISmoled l|2009h . 
With these convective parameters, overall properties of 



Table 2. Convective parameters considered for non-modulated 
initial models. Parameters a B , o c , a^, a p and 7 r are given in the 
units of standard values (o s = a c = 1/2-^2/3, a j = 8/3\/2/3, 
a p = 2/3 and 7 r = 2y/3; see lSmolec fc Moskalikl l2008al. for de- 
tails). 

a » m » s a c a d a p a t 7r 
1.5 0.65 1.0 1.0 1.0 0.0 0.0 1.0 



the Galactic Cepheids, both fundamental mode and first 
overtone pulsators, were modelled successfully. Here, we 
only increased the eddy- viscous dissipation (<a m parame- 
ter) in order to match the typical pulsation amplitudes 
of the RR Lyrae stars. We note that the modelling of 
RR Lyrae light curves is a difficult problem and some 
discrepancies with the observati ons still remain (see e.g. , 
th e detailed comparis on done by iKovacs fc Kanbuij (Il9981 ) 
or iFeuchtingej (|l999h V Nevertheless, the general proper- 
ties of the RR Lyrae light curves are quite well repro- 
duced with the current convective hydrocodes. In particular, 
the light curves of individual stars can be nicely matched 
with the hydrodynamic m odels (see e.g., iMarconil |2009| ; 
iMarconi fc Clementinill20051 ). 

The static models, the first step in our modelling, con- 
sist of 150 mass zones extending down to 2-10 6 K. Fifty outer 
zones have equal mass, down to the anchor zone in which 
the temperature is set to 11000 K (hydrogen ionisation). The 
mass of the remaining zones increases geometrically inward. 
In the second step of our modelling procedure, we have com- 
puted non-linear full amplitude models. To reach the limit- 
cycle pulsation, we have computed 2000 pulsation cycles by 
default. These non-linear non- modulated models were sub- 
sequently used as initial models for model integrations with 
variable turbulent convection. 

Several sequences of modulated models were computed. 
The properties of these models are summarised in Table [3] 
We investigated the effects of different amplitudes of the 
turbulent convection modulation, different modulation pe- 
riods, and different modulation shapes (see Fig. [I}. In each 
set, several models of different effective temperatures, ex- 
tending through the whole instability strip, were computed. 
For most of the computed models, the physical parameters 
of set A were used (Table[T|. We note that the mean physical 
parameters of the modulated models, such as mean radius or 
mean effective temperature, are affected by the modulation 
of turbulent convection and differ from the static equilibrium 
values. In the following, we will refer to the particular mod- 
els by providing the adequate set name from Tableland the 
effective temperature of the underlying static model, T| ff . 



3 RESULTS 

In this section we analyse the computed model sequences, 
described in Section 12.31 focusing on the light curve vari- 
ation through the Blazhko cycle, and on the properties of 
the frequency spectra. First, in Section [3]T] we describe the 
overall properties of the light curve modula tion in terms 
of th e Fourier decomposition parameters fe.g.. lSimon fc Led 
Il98ll) . The variety of behaviours and trends we compute 
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Table 3. Parameters of turbulent convection modulation for the 
model sequences considered in this study. For definitions of A, te 
and tf see Fig. [T] 



Set 


physical 


function 


A 


tn 


t{ 




params . 








Ml 


A 


sine 


10% 


60 d 




M2 


A 


sine 


20% 


60 d 




M3 


A 


sine 


30% 


60 d 




Ml 


A 


sine 


20% 


40 d 




MS 


A 


sine 


40% 


60 d 




M6 


A 


sine 


50% 


60 d 




M7 


A 


sine 


50% 


40 d 




M8 


A 


sine 


50% 


120 d 




M7L4 


B 


sine 


50% 


40 d 




M7L6 


C 


sine 


50% 


40 d 




M7L7 


D 


sine 


50% 


40 d 




MT1 


A 


triang. 


50% 


60 d 


45 d 


MT2 


A 


triang. 


50% 


60 d 


15 d 



can serve for comparison with continuous observations of 
Blazhko variables by satellite missions. In Section 13.21 we 
compare our results with Kepler observations of RR Lyr, 
for which a detailed analysis of t he light variation was al- 
ready done (|Kolenberg et al]|201ll ). In Section[3]3]we briefly 
analyse the pulsation period changes in our models and in 
Section f3. 41 we present the results of the frequency analysis 
of our models, focusing on the overall properties of the fre- 
quency spectra and how these compare with properties of 
Blazhko variables and RR Lyr in particular. Finally, in Sec- 
tion [33] we provide a more detailed discussion on the period 
doubling phenomenon. 

3.1 Light curve variation through the Blazhko 
cycle 

To study the light curve variation through the Blazhko cycle 
we u se the time dependent Fourier analysis (|Kovacs et al.l 
1987). We fit the Fourier series to the light variation of the 
consecutive pulsation cycles, 

m — A + 2^ Aj sin(jojt + tpj), (1) 

3 

and analyse time variation of the Fourier decomposition pa- 
rameters of low order, the amplitude, Ai, amplitude ratio, 
R21, and phase difference, <p2i, 

Rz\=Ai/Ai, tp2i = <p2-2(p!. (2) 

We note that the derivation of instantaneous amplitudes and 
phases through the time-dependent Fourier analysis is jus- 
tified, as the modulation is slow compared with the period 
of oscillation. 

For the analysed smooth hydrodynamical data, the 
Fourier parameters also vary smoothly with time. If period 
doubling occurs, the alternating shapes of the consecutive 
pulsation cycles manifest in series of wiggles that appear in 
the run of the Fourier parameters. In Figs. (3}{5]we plot dif- 
ferent relations between the Fourier parameters for different 
models. We present the effects of different amplitudes of the 
mixing-length modulation (Fig. [3}, different modulation pe- 



riods (Fig. [4} and different modulation shapes (Fig. [5} on 
the light curve variation. The instantaneous mixing length, 
a, R21 and ^21, are plotted versus the amplitude, Ai, in the 
top, middle and bottom rows of these figures, respectively. 
As in our models consecutive Blazhko cycles are repetitive, 
relations take the form of closed trajectories. In order to vi- 
sualise the variation across the instability strip, we display 
the models with different effective temperatures, for which 
we choose T c s ff =6300K (left columns), T c s ff =6600K (mid- 
dle columns; models to be compared later with RR Lyr) 
and TJfj =6800 K (right columns; models with period dou- 
bling). All models are characterised by the physical param- 
eters of set A (Table [!}. On each trajectory in these figures, 
a plus sign indicates the phase of maximum mixing length, 
and the direction of the time flow along the trajectory is 
depicted schematically. In addition, in each panel crosses 
connected with dashed line indicate the location of several 
non- modulated models with fixed mixing-length values, in a 
range covered by our modulated models. Below we describe 
the properties of the light curve variation connected with the 
location within the instability strip and different parameters 
of the modulation of the turbulent convection. 

3.1.1 General properties of the light curve modulation. 
Variation through the instability strip. 

Analysis of Figs. [3}[5] reveals some general properties of 
the computed trajectories, independent of the parameters 
of turbulent convection modulation (amplitude, period and 
shape). A basic observation is that the higher the effective 
temperature is, T| ff , the higher is the mean pulsation am- 
plitude, A\. Also, the higher the mean pulsation amplitude, 
the smaller the range of variation of R21 and 9521 ■ For the 
models with 7^=6300 K we are close to the linear red edge 
of the fundamental mode. In fact, for the models with a 
strong modulation of the turbulent convection (e.g., of the 
sets M3 and M6 in Fig. |3J), phases of large instantaneous 
mixing-length values (a > 1.65) correspond to the pulsa- 
tionally stable equilibrium models. On the other hand, at 
higher temperatures (7^=6800 K), small values of the mix- 
ing length (a ^ 1.05) at some modulation phases, if adopted 
in non-modulated models, would lead to first overtone (IO) 
pulsation (see the long-dashed lines in the figures). It is 
clear that the temporary light curve of the model in which 
turbulent convection is modulated is significantly different 
from that of the non-modulated model adopting the same 
value of the mixing length. In addition, hysteresis is clearly 
visible for the modulated models. Very different shapes of 
the light curve are possible at two phases characterised by 
the same instantaneous mixing length. The described fea- 
tures are in qualitative agree ment with the observations. 
lJurcsik. Benko &; Szei d3 i|2002l ) noted that the light curves 
of Blazhko stars are never (at any phase of the Blazhko cy- 
cle) like those of non-Blazhko stars, thus always distorted. 
A detailed an alysis of the light curv e variation in Kepler 
RR Lyr data (|Kolenberg et al.ll201ll ) reveals a similar be- 
haviour of the Fourier parameters as plotted in the figures. 
A more quantitative comparison is postponed to the next 
section. 

Another interesting feature is visible in the a vs. Ai 
plots (top rows of Figs. [3H5}. We note that the minimum 
value of the amplitude does not coincide with the maximum 
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Figure 3. Light curve variation through the Blazhko cycle for the models with different amplitude of mixing-length modulation, 10 
per cent (set Ml), 30 per cent (set M3) and 50 per cent (set M6). The instantaneous mixing-length, R21 and 9221 are plotted versus 
Ai in the top, middle and bottom rows of the Figure. In consecutive columns we plot the results for models of different T B „ (6300 K, 
6600 K and 6800 K). A plus sign on each trajectory indicates the phase of maximum mixing length. The direction of the time flow along 
each trajectory is shown on a schematic circle. Crosses connected with a long-dashed line correspond to the non-modulated models with 
different values of the mixing length. For Tf ff =6300K and a > 1.65 the fundamental mode is linearly stable and the models were not 
computed. For 7^=6800 K and a ^ 1.05 the models switch into first overtone (IO) pulsation (two models in the Figure). 



value of the mixing-length parameter, as one might naively 
expect, but it is delayed. The reason for such a delay is not 
clear. It shows the inertia of the dynamical system, which 
does not adjust instantaneously to the changing dissipation. 

Considering the relations R21 vs. Ai and <^2i vs. Ai 
(middle and bottom rows of Figs. [3] [5j) we see that for low 
temperature models (of small mean pulsation amplitudes) 
trajectories have the shape of double loops, while at higher 
temperatures (and high mean pulsation amplitudes) single 
loops are present. The direction of trajectories (in case of 
double loops, the direction of the larger loop) is always 
clockwise for 1^21 vs. Ai trajectory (bottom rows of Figs. [3]- 
[5]). For R21 vs. Ai trajectories (middle rows of Figs. [3]^5j) 
the direction is either counter-clockwise (for cooler models, 
T c s ff =6300K and T e s ff =6600 K) or clockwise (for the hottest 
models, =6800 K). The analysis of additional models 
across the instability strip reveals that the transition occurs 
at around T^ s =6700 K. For these models we deal with an al- 
most single- valued dependence of R21 on Ai . We note that 
the direction of the ip2i vs. Ai loop is related to a partic- 
ular asy mmetry of the triplet components in the frequency 
spectra (|Szeidl fc Ju rcsik 200fJ). We discuss this problem in 
more detail in Section [3.4. 21 

Considering the 7?2i vs. Ai relation, we first note that 
intuitively we expect that R21 should correlate with Ai. The 
higher the amplitude, Ax, the larger the non-linearity and 



consequently larger the contribution of the harmonic terms 
(and thus, higher i?2i). In our models, in which both R21 
and Ai are functions of time, the relation between R21 and 
Ai is complicated. For the hottest models, of higher mean 
pulsation amplitude, R21 correlates with Ai. For the lower 
temperature models, of lower mean pulsation amplitudes, 
the trajectories are more extended ("blown-up") and the 
relation is not obvious. At some phases in the Blazhko cycle, 
a high amplitude Ai is accompanied by lower values of R21 ■ 

Period doubling is clearly visible for the highest tem- 
perature models (T* s = 6800 K) displayed in Figs. [3HSI m 
which the mixing-length modulation is strong (Fig. [3]). Dur- 
ing the phases with period doubling, the trajectories are no 
longer smooth. The series of wiggles is present, which indi- 
cates that the consecutive pulsation cycles alternate. Period 
doubling will be discussed in more detail in Section f3. 51 

In addition to the models running horizontally in the 
instability strip, we have computed some models with dif- 
ferent M/L ratio. Four sequences of models were computed, 
M7L4, M7, M7L6 and M7L7 (see Table |3j), in which the lu- 
minosities vary from 40 L© to 70 Lq. The mass is the same 
in all these model sequences (0.65 Mq) and T| ff was varied. 
In all these sequences the amplitude of the mixing-length 
modulation was 50 per cent, the modulation period 40 days, 
and the shape of the modulation was sinusoidal. Except for 
the shifts in the computed trajectories, connected with the 
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Figure 5. The same as Fig. [3] but for the models with different shapes of mixing-length modulation; sinusoidal (set M6) and asymmetric 
triangular (MT1 and MT2). 
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different luminosities of the initial models, the trajectories 
are very similar and hence not presented in a separate figure. 
We have noted that the direction of the R21 vs. A\ trajecto- 
ries for high temperature models depends on the luminosity 
and is clockwise for the low-luminosity models (40 L0 and 
50 Lg) and counter-clockwise for the high-luminosity models 
(60 L© and 70 Lq), in which R21 displays a flat dependence 
on Ai . Consequently, a clockwise direction of the R21 vs. Ai 
trajectory is present only in higher-temperature and lower- 
luminosity models. 

3.1.2 Different strength of the turbulent convection 
modulation. 

In Fig. [3]we plot the results for the models with a different 
amplitude of the turbulent convection modulation, A — 10 
per cent (Ml), A = 30 per cent (M3) and A = 50 per 
cent (M6). The period and shape of the modulation are the 
same for these three sets (Pb = 60 d, sine). As expected, 
the stronger the modulation, the more complicated ('non- 
linear') trajectories we get and higher the ranges of variation 
of the Fourier parameters. This property of the computed 
trajectories will be used in the next section to estimate the 
required amplitude of the mixing-length modulation to re- 
produce the RR Lyr observations. We also note that, de- 
pending on the strength of the modulation, the mean values 
of Fourier parameters vary. This concerns the amplitude, Ai, 
and the amplitude ratio, R21. The mean phase difference is 
not very sensitive to the strength of the modulation. For 
the most convective models of T c s ff =6300K, the stronger the 
modulation, the smaller the amplitude and the amplitude 
ratio. 

3.1.3 Different period of the turbulent convection 
modulation. 

In Fig. [4] we plot the results for models with different pe- 
riod of turbulent convection modulation, Pb = 40 d (M7), 
P B = 60 d (M6) and P B = 120 d (M8). The amplitude and 
shape of the modulation are the same for these three sets 
(A — 50 per cent, sine). The basic observation is that for 
longer modulation period, the computed loops are larger 
(more "blown- up"), and thus, the larger the ranges of varia- 
tion of Fourier parameters, but the effect is not as strong as 
in case of different amplitudes of turbulent convection modu- 
lation. We also note that the range of effective temperatures 
in which the period doubling occurs, depends on the mod- 
ulation period, but we postpone the detailed discussion to 
Section 13.51 

3.1.4 Different shape of the turbulent convection 
modulation. 

In Fig. [5] we plot the results for models of sets M6 (sinusoidal 
modulation of the mixing length), MT1 and MT2 (asym- 
metric triangular modulation). The period and amplitude 
of mixing-length modulation are the same for these three 
sets (Pb = 60 d, A — 50 per cent). We conclude that the dif- 
ferences between the trajectories are rather minute and the 
light curve variation is very similar for the different shapes 
of modulation. 



3.2 Comparison with RR Lyr 

RR Lyr is not only the prototype of the RR Lyrae vari- 
ables but also a famous Blazhko variable with a strongly 
modulated light curve. Its Blazhko period is subject to 
small variations, its p resent value being 39.1 ± 0.3 days 
(jKolenberg et al.ll2011f ). Thanks to its location in the Ke- 
pler field of view, we now have excellent, nearly continuous 
photometric data covering slightly more than three Blazhko 
periods - seasons Q1 +Q2 of the Kepler long cadence data 
|jenkins et alj|2010al lbh of which Ql is already publi(0. 

In this section we compare the light curve variation 
through the Blazhko cycle, as it is observed in RR Lyr, 
with our models. The Kepler photometric passband is rather 
wide, cover ing the combined range of the standard V and R 
passbands (|Koch et al.ll2010T ). It justifies the direct compar- 
ison of our model bolometric light curves with the Kepler 
RR Lyr light curve. This simplification does not affect the es- 
timate s presented in this section. We note that lNemec et ail 
l|201ll ) derived the transformation between the Fourier pa- 
rameters in Kp and in V passbands based on the photometry 
of three non-modulated RR Lyrae stars in both bands. The 
differences are systematic but very small. 

In Fig. [6] we plot the 1^21 vs. A\ and P21 vs. A\ rela- 
tions as observed for RR Lyr (season Q2 of Kepler data). 
The ranges of variation of the Fourier parameters are rather 
large. In our models, as noted in the previous section, the 
larger the amplitude of turbulent convection modulation is, 
the larger is the range of variation of the Fourier parameters. 
Now, we can estimate the strength of mixing-length modula- 
tion necessary to explain the RR Lyr observations, through 
comparing the model and observed ranges of variation of 
^4i, R21 and f>2i- To this aim, for each Fourier parameter c, 
c € {Ai, R21, <p2i}, we define the parameter Ac, 

a C max ~~ C m in , v 

Ac = 2 , (3) 

Cmax ~T Cmin 

constructed using the minimum, c m i n , and maximum, c max , 
values of c during the modulation cycle. We note that the 
value of Ac is independent of simple scaling of the Kepler 
data by a constant factor. For RR Lyr we have AA^ = 
0.62, AP21 = 0.38 and A</? 2 i = 0.27 jKolenberg et al.ll201ll . 
Fig. O . The computed values for the models with a different 
strength of the turbulent convection modulation, A = 10 per 
cent (set Ml), A = 30 per cent (set M3), and A = 50 per 
cent (sets M6 and M7) are collected in TableUJ For each set, 
the values for seven models of different are computed. 
The values that agree within 20 per cent with RR Lyr values 
are marked with an asterisk. Note that the sets M6 and M7 
have the same amplitude of the mixing-length modulation 
(A = 50 per cent), but for set M7 the modulation period is 
shorter (40 days, very close to RR Lyr's modulation period). 

To reproduce the ranges of the Fourier parameter vari- 
ation in RR Lyr, a mixing-length modulation with an am- 
plitude equal to at least 30 per cent is necessary. The best 
match is found for an amplitude of the mixing-length mod- 
ulation equal to 50 per cent (sets M6 and M7). Then, 
both Aj4i and AP21 can be matched for the model with 
T c s ff =6600K. For models with T c s ff =6600K, the mean pulsa- 
tion period (w 0.573 d) matches RR Lyr's period (0.567 d, 

3 http:/ /archive. stsci.edu/kepler/ 
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Table 4. Ranges of variation of the Fourier parameters for the different models of sets Ml, M3, M6 and M7 compared with the RR Lyr 
Kepler data. In the first column effective temperature of the initial non- modulated model, is given. In the consecutive columns the 
ranges of variation of Fourier parameters, AAi, AR21 and Aip2l (see eq. [3j are given for different sets of modulated models. In the last 
row, the data for RR Lyr are provided for reference. Model values that agree with the corresponding RR Lyr value within 20 per cent 
are marked with an asterisk. 
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Figure 6. Fourier parameters, tp2i and R21, plotted vs. Ai for 
Kepler RR Lyr data (season Q2) and for two models of set M3. 



iKolenberg etTal] l201ll ) almost exactly. The range of the 
phase variation, Ay>2i, however, is slightly higher for these 
models than is observed in RR Lyr. 

The most stringent constraints on the Stothers model 
come from highly modulated stars, like RR Lyr. Tiny mod- 
ulations as observed for many Blazhko stars can be easily 
reproduced assuming a small amplitude of turbulent con- 
vection modulation in our models. It is now evident that, 
if the mechanism proposed by Stothers is responsible for 
the Blazhko effect, the observed large modulations of the 
light curves (as we observe in RR Lyr) require significant 
changes of the mixing-length over the Blazhko cycle. Note 
that the modulation amplitude, A, as defined in Fig.[TJ is ac- 
tually a semi-amplitude. For the considered models, A — 50 
per cent means that the mixing length, a, varies in a huge 



range, from 0.75 to 2.25. In our opinion, such large changes 
in the effectiveness of the turbulent convection on relatively 
short time-scales of typically tens to hundreds of days (typ- 
ical Blazhko periods) are highly unlikely (see discussion in 
Section [5]). 

Now we focus on the shape of the Fourier parameter 
trajectories. A close inspection of Figs.|3}[5j particularly the 
models with 1^=6600 K (for which the period of the fun- 
damental mode matches that of RR Lyr) shows that the 
model Lp2i vs. Ai loops are quite similar to the RR Lyr loop 
(top panel of Fig. [SJ). What we typically see in our models 
(around T| ff =6600 K) is a larger ("blown-up") part of the 
loops on the left side and cusps on the right side. This is 
what we observe in RR Lyr, however, the cusp in our mod- 
els is located at higher values of if 21 compared to RR Lyr. 
A comparison of Kepler RR Lyr data with two models of set 
M3 (for which we get the best Ai^2i match) is shown in the 
upper panel of Fig. [6] We note that also the direction of the 
time flow is the same for both RR Lyr and for the models 
(clockwise). 

Considering the i?2i vs. Ai loops (bottom panel of 
Fig. |SJ) we cannot reproduce the behaviour that we observe 
in RR Lyr. In RR Lyr 7?2i anticorrelates with Ai through 
the majority of the Blazhko cycle. This is not reproduced by 
our models, although at some phases of the modulation R21 
increases as Ai decreases. The bottom panel of Fig. [6] pro- 
vides a direct comparison of the model (set M3) and RR Lyr 
loops. Also here, the direction of the time flow is the same for 
our models and for RR Lyr (bigger loop, counter-clockwise). 



3.3 Changes of the pulsation period 

In this section we analyse the changes of the pulsation pe- 
riod. In the mathematical description, a period change is 
equivalent to a change of the pulsation phase. From the ob- 
servational point of view, the two are indistinguishable. It 
is the physical interpretation where the difference occurs. 
Period changes are caused by the overall changes of the stel- 
lar structure, while phase changes may result, e.g., from the 
nonlinear interaction of pulsation modes (which does not 
affect the stellar structure). 

In our models, due to the changes of the convective 
structure, the pulsation period changes during the modula- 
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Table 5. Amplitudes of the period change, SP/P, for several 
models of different T s „ and different amplitudes of turbulent con- 
vection modulation, 10 per cent (Ml), 30 per cent (set M3) and 
50 per cent (set M6). 
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tion. The amplitude of the period variation, SP/P, is most 
sensitive to the amplitude of turbulent convection modula- 
tion, and hence can be used to estimate the required strength 
of mixing-length modulation in a similar way we have done 
in the previous section. In Table[5]we list the period changes 
for the models with different amplitudes of the turbulent 
convection modulation, A = 10 per cent (set Ml), A — 30 
per cent (set M3) and A = 50 per cent (set M6). The pe- 
riod changes were estimated based on the radius variation 
of our hydrodynamic models, using the time difference be- 
tween consecutive radius maxima as an estimate of the in- 
stantaneous period. Then, the full amplitude of the period 
variation was used to derive SP/P values. It is evident that 
the larger the amplitude of mixing-length modulation is, the 
larger the amplitude of period variation. Also, the trend of 
decreasing amplitude of the period variation with increasing 
effective temperature is clear. In Blazhko RR Lyrae vari- 
ables, the pulsation period changes are clearly detected, with 
a typical amplitude, SP/ P, between 0.2 and 1.4 per cent 
l|Molnar fc Kollathl I2010TI. For RR Lyr the period change 
is 0.83 per cent ( Kolenberg et all 1201 lj ) . Comparison with 
values collected in Table [5] indicates that in order to repro- 
duce such period variations, large amplitudes of the mixing- 
length modulation are necessary, at least of the order of 50 
per cent, in agreement with the value deri ved in the previous 
sectio n and in agreement with estimate of lMolnar fc Kollathl 
(|2010l V based on the analysis of the periods of linear equi- 
librium models. 



3.4 Analysis of frequency spectra 

In this section we focus on the frequency anal ysis of our mod- 
els, fo r which we use the Period04 package |Lenz fc Bregerl 
|2005| ). Differently from SectionPXTl we now analyse the light 
variation over many pulsation periods and many Blazhko cy- 
cles. At this long time-scale the system is stationary and we 
can fit the data with a Fourier sum of the following form: 



N , 

A + ^ \ A k sin (27rfe/ i + <pk) + 

k = l ^ 

[27r(fc/ + / B )t + 4 
2w(kf -f B )t + 4>k 



AJ. sin 



A k sin 



A{ + sin 2tt (kfo + ifa)t + <j>\ 

J 

Air sin 2n(kfo-ih)t + 4 

L 

Bj sin (2njf B t + 4>bj) ■ 



+ 



+ 



(4) 



3=1 



In the above formula, the first line corresponds to the main 
pulsation frequency, fo, and its harmonics, kfo, the second 
and third lines correspond to the components of the triplets, 
fc/o±/B, the next two lines describe other higher-order mul- 
tiplet components and in the last line we account for the 
modulation frequency, /b, and its harmonics, jf B - We as- 
sume that the components of the multiplets are equidistant. 
Also the modulation frequency, /b, is known, as it is equal 
to the inverse of the assumed period of mixing-length mod- 
ulation (see Table [3]). Consequently, only one frequency is 
determined from the data, the main pulsation frequency, 
which is not known a priori (the non-linear period differs 
from the linear one). By default we use N = 19, J = 5 and 
L = 2 for all the models (note that signal at 2/b is strong 
in our models). 

The length of the hydrodynamic data we analyse corre- 
sponds to four full Blazhko cycles. To speed up the compu- 
tations, the sampling of the model light curve is degraded to 
roughly 60 points per pulsation period (~twice the Kepler 
resolution for RR Lyr). 

Below, we present some details of our analysis for one 
particular model, and later (Section I3.4.2p we discuss the 
general properties of the frequency spectra of our models 
and compare them with the observations. 



3.4-1 Frequency spectra analysis - particular case 

We analyse one particular model from set M6 with T| ff = 
6800 K. The model displays a clear period doubling in the 
computed light curve (see, e.g., Fig. [3j . The mean pulsation 
period for this model is Pi = 0.517 d, which is not far from 
the RR Lyr period (Pi = 0.567 d, Kolenberg et al]|201ll ). 
As noted before, the best match with RR Lyr's period is 
obtained for the model with T < f ff =6600K, but this model 
does not display the period doubling. 

For the discussed model, the variation of the bolometric 
magnitude with time for slightly more than one Blazhko 
cycle (70 days) is shown in the upper panel of Fig. [5] Period 
doubling is clearly visible during the phases of maximum 
pulsation amplitude. The lower panel of Fig. [2] shows these 
phases. The period doubling is obvious and lasts for several 
pulsation cycles. However, its amplitude is not very large. 

In Fig. [7] we plot some results of the frequency anal- 
ysis. The upper panel of Fig. [7] shows the close-up around 
the fundamental mode frequency, fo, and its harmonic, 2/o, 
after prewhitening the data with the fundamental mode fre- 
quency and its harmonics (dashed lines), and five consecu- 
tive components of the multiplet structure. In the middle 
panel of Fig. [7] we show the vicinity of the fundamental 
mode frequency. All removed frequencies are marked with 
dashed lines. Clearly, many more higher-order side peaks 
are visible. Residual power at the position of fo corresponds 
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Figure 7. Frequency spectrum for the particular model of set 
M6 (T*~ = 6800 K) after prewhitening with the Fourier sum de- 
scribed by Eq. JJJ (upper panel). Two close-ups are shown: at the 
location of fundamental mode frequency (middle panel) and at 
its subharmonic, 1/2 fo (bottom panel). Removed frequencies are 
marked with dashed lines. 



to the long-term evolution of the model toward the final 
limiting cycle pulsation (which is of exponential character). 
The bottom panel of Fig. [7] shows the vicinity of the half- 
integer frequency (HIF), 1/2 fo- The sign al at this frequenc y 
is a signature of the period doubling (see lSzabo et alJl2010h . 
which appears in the discussed model during a short fraction 
of the modulation cycle. Consequently, the signal is strongly 
modulated with the Blazhko period and hence many equally- 
spaced peaks are clearly visible, with the separation corre- 
sponding to the modulation frequency. The envelope of the 
signal located at the HIF's is very regular and resembles a 
Gaussian. The width of this Gaussian corresponds to the 
duration of the period doubling behaviour observed in the 
model (6 — 7 days). 

Our analysis is focused on the first-order side peaks, the 
triplets. The relations between their amplitudes, AT" , A^ and 
A+ (see Eq. 0} were studied for several Blazhko stars ob- 
served from the ground (e.g.. Ijurcsik et al.ll2006h as well a s 
for the Kepler RR Lyr observations jKolenberg et alj|201ll ). 
For our models, the relations between the amplitudes of the 
signals at different frequ encies are shown in Fie . [8] This is 
the analog of fig. 7 from iKolenberg et alj (|201lf ). where the 
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Figure 8. Amplitude ratios, Afc/Ai, AZ /A~ , A+ /A+ and ampli- 
tude of the half-integer frequencies, plotted versus the frequency 
(harmonic orders indicated by vertical dashed lines). Results for 



the model of set M6 with T s , 



6800 K. 



results for RR Lyr are presented. The model we discuss now 
has physical parameters and a mean period close to RR Lyr, 
but the modulation period is slightly longer (60 days) . How- 
ever, the described features of the frequency spectrum do not 
depend strongly on the modulation period and the analysis 
below is comparative. 

For Ak/Ai we observe an exponential decrease with 
harmonic order k. At k = 7 the amplitude ratio drops to 
very small values (just like in RR Lyr). For the modulation 
components the decrease of amplitude ratios (AT" /AT and 
A^ /Af) is less steep. It is also less steep for A\Z /AT/ than 
for A~/"/A+ - all in agreement with RR Lyr. For A\Z jA~ we 
observe a tail at high k also in agreement with RR Lyr. For 
low k, At; /A 7 increases (above 1) in the case of RR Lyr. 
Here it is not the case, but such behaviour can be obtained 
for other models (see next section). 

For the amplitudes of the HIF's resulting from period 
doubling, in the discussed model the highest peak is located 
at 1/2 f followed by 3/2 f and 5/2 /q. For RR Lyr the high- 
est peak is observed at 3/2 fo followed by 5/2 fo and 1/2 /o. 
We discuss this discrepancy in more detail in Section \3. 51 



3.4-2 Frequency spectra analysis - model sequences 

A detailed frequency analysis was conducted for all the 
model sequences discussed in this paper. Typical results for 
models with T , < f ff =6800 K and different parameters of turbu- 
lent convection modulation are shown in Fig. [5] These mod- 
els are similar to that of set M6 discussed in the previous 
section, have the same amplitude of mixing-length modula- 
tion (50 per cent), but a longer modulation period (set M8) 
or a different shape of the modulation (triangular for MT1 
and MT2). We stress that the results discussed below are 
characteristic for all our model sequences. 

For the amplitudes of the harmonic frequencies, A^ /Ai 
(panel a of Fig. [9} , an interesting behaviour is observed at 
harmonic order around 7. The amplitudes are already very 
small at this order so the plots are in logarithmic scale to 
reveal the details. Local minima are clearly visible. They 
fall at k = 7, except for the model with a longer modulation 
period (set M8, k = 9). In all these models period doubling 
is present, however the origin of the discussed minima and 
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Figure 9. Relations between the amplitudes of the triplet com- 
ponents for the models of sets M6, M8, MT1 and MT2 (all with 
1^=6800 K). Panel (a): A^/A\ vs. k in logarithmic scale; panel 
(6): A^/A^ vs. k; panel (c): A^/A^ vs. k, and panel (d): Qk 
vs. k. The distribution of Q values from the MACHO database 
peaks at 0.3, which is marked with the horizontal dashed line for 
reference. 



their possible connection with period doubling is not clear 
(see also Section I33)) . 

Now we discuss the amplitudes of the modulation com- 
ponents, A^/Ai and A^/Af (panels b and c of Fig. 0). 
For A^ /A~[ we observe more 'erratic' behaviour than for 
A^/Af. At low orders, an increase of the amplitude ratio 
(above 1) is possible, just as it is observed in RR Lyr. A~^ j A\ 
decreases with increasing order. A plateau is visible for har- 
monic orders between 3 and 5. 

Panel d of Fig. [§] shows the va riation of the Qk param- 
eter defined as (jAlcock et al.ll2003r i. 



Qk 



■a: 



At + A: 



(5) 



with k. The line at Q = 0.3 is marked for reference (the 
distri bution of Q from the MACHO database p eaks at this 
value, lAlcock et al. l l2003l ; lKolenberg et alj|201lh . It is clear 
that in all our models Qk is positive at low harmonic orders 
and hence, the amplitudes of the higher frequency triplet 
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Figure 10. Qi plotted vs. the Tf ff for models with different lu- 
minosities of the initial non-modulated models (sets M7L4, M7, 
M7L6, M7L7). 



components are higher. In the intermediate range of har- 
monic orders the situation is reverted and at largest k the 
higher frequency components are higher again. 

We focus our attention on the triplet components 
around the main frequency fo. The amplitudes of these 
triplet components are the most robust outcome of our 
model and should be compared with observation first. In 
about 75 per cent per cent of Blazhko RR Lyrae stars, 
the higher frequency side peak has a larger amplitude than 
the lower frequenc y side peak and thus Q\ is positive 
l|Alcock et al.1120031 '). This is in agreement with our models. 
However, in 25 per cent of the observed Blazhko variables 
Qi is negative, as the lower frequency side peak has a larger 
amplitude. Unfortunately, in none of our models this is the 
case. In Fig. \W\ we plot the values of Qi vs. T| ff for the 
models of sets M7L4, M7, M7L6, M7L7, with different lu- 
minosities of the initial non-modulated models. The trend 
of the decreasing Q\ values with decreasing luminosity as 
well as with decreasing temperature is clearly visible. How- 
ever, in no case Q\ is negative. This is true for all model 
sequences we have investigated, regardless of the properties 
of turbulent convection modulation (period, amplitude and 
shape). As in a quarter of known Blazhko RR Lyrae stars 
Qi is negative, we regard the failure to reproduce negative 
Qi values in our models as another significant challenge for 
the Stothers mechanism. 



3.5 Period doubling phenomenon 

The period doubling phenomenon, which manifests in alter- 
nating shapes of the light curves of the consecutive pulsa- 
tion c ycles, was discovered v ery recently in Kepler RR Lyr 
data ( Kolenberg et al]|2010al ). The effect is also clearly vis- 
ible in two other Kepler targets, V808 Cyg (KIC 4484128) 
and V355 Lyr (KI C 7505345), both being Blazhko variables 
l|Szab6 et alj|201ol ). In some other Blazhko stars, the confir- 
mation of period doubling requires longer observation. The 
strength of period doubling depends on the phase in the 
Blazhko cycle. Period doubling was not found in any non- 
modulated RR Lyrae variable so far (|Szab6 et alj|20ich . 

The phenomenon of period doubling is not new. Period 
doubling is clearly present in RV Tauri stars which show 
alternating deep and shallow minima in their light and ra- 
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dial velocity curves. It was also found in radiative hydro- 
dynamic models of W Vir stars |Buchler fc Kovacsl 1 19871 ; 
iKovacs fc Buchlerl ll988T), Cepheids and BL Her variable s 
jMoskalik fc Buchlerlli990l,ll99ll ; lBucIiler fc Moskalik|[l993) . 
iMoskalik fc Buchlerl (jl990l ) traced the origin of the period 
doubling in hydrodynamic models to the destabilising role 
of the half integer resonance (2n + l)uo = 2uik between the 
fundamental mode and a higher order overtone (n equal to 1 
or 2). They showed that the period doubling can occur close 
to the resonance centre when the reson ant overtone mode is 
either excited or only weakly damped. Mos kalik fc Buchlej 
(1990) also showed that, depending on the amount of dissi- 
pation in the model, the half-integer resonance leads either 
to single period-doubling or to the period-doubling cascade 
(Feigenbaum cascade) and chaos. 

The period doubling phenomenon in Blazhko RR Lyrae 
stars observed b y Kepler was analysed in detail by 
ISzabo et al.l l|2010l ). who also proposed the underlying mech- 
anism. Their convective hydrodynamic models (with fixed 
convective parameters) display period doubling, the origin 
of which was traced to the 9:2 resonance between the fun- 
damental mode, and the ninth order overtone, 9ujo — 2lj$. 
In these mode ls, the ninth overtone is a trapped envelope 
mode (see e.g jBuchler fc KollathlkoOll ) with a much higher 
growth rate than that computed for the neighbouring over- 
tone modes. As noted above, the weak damping of the ninth 
overtone favours the occurrence of per iod doubling. In their 
detailed a nalysis Kollath et al.l l|201 ll ) used the relaxation 
technique ('St ellingwerilll974l ) to determine the stability of 
the fundamental mode pulsation through the Floquet stabil- 
ity coefficients. They showed that indeed the fundamental 
mode is destabilized through the 9u>o = 2ug resonance and 
excluded all other possible half-integer resonances. In addi- 
tion, they showed that the 9c^o = 2u>g resonance can lead not 
only to a single period doubling, but also to the period-four 
and period-eight limit cycles (the Feigenbaum cascade). 

Our models also display period doubling, clearly visi- 
ble in the light curves (Fig. [2] Figs. [3}{5]l and in the fre- 
quency spectra (Figs. [7]and[8]|. It occurs only in models 
of higher temperatures, with T c s ft in a range 6700 K-6900 K 
(depending on the model sequence) and always at phases 
around maximum pulsation amplitude and minimum values 
of the mixing length (see right colum ns of Figs-ErEp- Ou r lin- 
ear analysis confirm s the findings of lSzabo et al.l ( 2010l ) and 
iKoIlath et alj (|201ll ). During the phases of minimum mix- 
ing length, in which period doubling occurs, the equilibrium 
models computed assuming corresponding, fixed values of 
the mixing length are very close to the 9ojo = 2u;g resonance 
centre. We note that the mixing length has to be sufficiently 
small. Period doubling does not occur in models with small 
amplitude of the mixing-length modulation (Fig. [3}. Also, 
the ninth overtone in our models is trapped in the outer 
envelope and is close to being unstable. 

For the period doubling to occur, the model has 
to be close to th e resonance condition. As discussed by 
ISzabo et al.l |2010l ) , if the Stothers mechanism is indeed op- 
erational in Blazhko variables and the convective structure 
of the star varies during the Blazhko cycle (as it does in our 
models), it is natural that period doubling occurs only at 
certain Blazhko phases, at which the physical conditions are 
favourable. Qualitatively, this is what we see in our models. 
Favourable conditions (proximity to the 9:2 resonance) oc- 
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Figure 11. Amplitudes of half integer frequencies for models of 
sets M6, M8, MTl and MT2 (initial non-modulated model was 
the same for all four models displayed, and its effective tempera- 
ture was 6800 K). 



cur during the phases of low mixing length. These are the 
phases at which period doubling occurs. We note that in 
our models period doubling is strictly repetitive and always 
appears at the same phase of the Blazhko cycle. In fact, at 
other phases the period doubling does not vanish entirely. 
Its remnants serve as a seed for the fast growth of alterna- 
tions during the next Blazhko cycle. We would like to stress 
that turbulent convection itself does not cause period dou- 
bling, which is a resonant phenomenon and may occur in 
purely radiative models as well. The modulation of turbu- 
lent convection in our models just changes the structure of 
the outer stellar layers, bringing the instantaneous model 
periods close to the resonance condition. 

As discussed above, the transient occurrence of the pe- 
riod doubling during the Blazhko cycle can be naturally in- 
terpreted within the Stothers model. Nevertheless, a more 
quantitative comparison with observations reveals some dis- 
crepancies. In our models period doubling is limited to the 
phases of maximum pulsation amplitude. For the three Ke- 
pler Blazhko RR Lyrae stars it is most prominent in the 
phases preceding the maximum pulsation amplitudes, but 
not only, it is also visible close to the phases of minimum pul- 
sation amplitude l|Szabo et al.ll201ol ). Certainly it is present 
in a much wider range of Blazhko phases than it is in our 
models. 

As for the frequency spectra, the amplitudes of the half 
integer components in all three stars observed by Kepler fol- 
low the same pattern. The highest peak is observed at 3/2 fo 
followed by 5/2 fo and 1/2 fo- This is not what we compute 
in our models, as already mentioned in Section 13.4.11 In 
Fig. [11] we show the amplitudes of HIF's, for models with 
different patterns of turbulent convection modulation (T^g 
is the same in these models, equal to 6800 K). These are 
the same models as displayed in Fig. [9] The highest peak is 
always located at 1/2 fo followed by 3/2 fo and 5/2 fo- For 
the latter two peaks the amplitudes are comparable. A local 
maximum is visible at 9/2 fo for our models of sets M6, MTl 
and MT2. Such a maximum is expected if the 9:2 resonance 
is indeed responsible for the period doubling. The fact that 
it is very weak and not present in all our models is a puzzle. 
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Figure 12. Light curve variation through the Blazhko cycle. R21 vs. Ai and ifi2i vs A\ for models in which different a-parameters of 
the convective model are modulated. 



4 ADDITIONAL MODELS 

The results presented so far reveal some disagreements be- 
tween the computed models and the observations. The fre- 
quency analysis shows that in all our models Qi is positive, 
so the higher frequency side peak around the fundamental 
mode frequency has a higher amplitude than the lower fre- 
quency side peak. In about 25 per cent of Blazhko variables 
Qi is negative. Also, we cannot reproduce the details of the 
light curve variation. In particular the anticorrelation be- 
tween amplitude ratio, R21, and amplitude, A\, observed in 
RR Lyr through the majority of the Blazhko cycle cannot 
be reproduced. Probably, the most important challenge for 
the Stothers mechanism is a huge range of variation of the 
mixing length. This large variation is required to reproduce 
the ranges of light curve variation in strongly modulated 
stars like RR Lyr, as well as the pulsation period changes 
we observe in Blazhko variables. 

To get an idea whether and how these difficulties may be 
overcome, we have computed additional sequences of models 
in which we vary other a-parameters of the turbulent con- 
vection model than the mixing length (see Section [2.11) . So 
far, we modulated the mixing length only, which is the natu- 
ral choice, as the mixing- length parameter controls the over- 
all strength of convection. Now, as an exercise, we vary the 
other parameters entering the turbulent convection model. 
We modulate the strength of eddy- viscous dissipation (a m ), 
the strength of the source function (a s ), the strength of 
the convective heat flux (q c ) and the strength of the tur- 
bulent dissipation (ay, turbulent-cascade term). During the 
mixing-length modulation, all the above terms were modu- 
lated simultaneously, as these a-parameters enter the model 
in pair with the mixing- length parameter (a a m , aa E , aa c 
and a^/a; see e.g.. lSmolec fc Mosk alik 2008a). We note that 
there is no physical justification behind the modulation of 
any particular term listed above. Only a detailed 3D magne- 
tohydrodynamical computation could clarify how a variable 
magnetic field could affect particular phenomena connected 



with the turbulent convection. Unfortunately, such compu- 
tations and even appropriate models do not exist currently. 

In all computed model sequences, the shape of modula- 
tion is sinusoidal, its period is 40 days and the amplitude of 
its modulation is 50 per cent. Particular a-parameters vary 
around the values defined in Table [2] 

In Fig. [T3] we plot the R21 vs. A\ and 9521 vs. A\ re- 
lations for models with three different initial temperatures 
(T c s ff =6300 K, 6600 K and 6800 K), just like in Figs.[3H5] The 
cross in each panel corresponds to the initial non-modulated 
model. Thick solid trajectories correspond to the model of 
set M7 in which the mixing length is modulated (default 
in this paper). As expected, the ranges of variation of the 
Fourier parameters are smaller if we vary other a-parameters 
of the model than the mixing length, as in each case only 
one term of the convective model is modulated. We focus 
on i?2i vs. Ai relation (upper panels of Fig. [12)) . Observed 
trends are not always simple (particularly for models with 
T e s ff =6300K) and they depend on the temperature of the 
model. The increase of R21 with decreasing Ai, as we ob- 
serve in RR Lyr, can be reproduced most easily when only 
the strength of the convective heat flux is modulated. Then, 
however, the range of variation of the Fourier parameters is 
small, despite a rather huge modulation of the convective 
heat flux. Therefore, it is difficult to propose a modulation 
to get a better model for RR Lyr. We have to vary more than 
one convective parameter to get a large range of variation 
of the Fourier parameters. For the model with =6600 K 
one could modulate only a s , a c and Od and keep the eddy- 
viscous dissipation fixed. On the other hand, Fig. [12] sug- 
gests that the same modulation adopted in the hotter model 
(7^=6800 K) would lead to increasing R21 with increasing 
Ai. 

We have also investigated whether the negative Q\ val- 
ues can be obtained through modulating particular terms in 
the convective model. In Fig. [13] we plot the Qi vs. T e s ff for 
the discussed models. Negative, albeit very close to zero, val- 
ues of Qi are obtained only for the hottest models for which 
either only convective heat flux is modulated (T° ff =6600 K 
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Figure 13. Q\ plotted vs. the T?„ for models, in which particular 
terms of the convective model were modulated. 

and 7^=6800 K) or only eddy-viscous dissipation is modu- 
lated (T c s ff =6800 K). Modulation of other components of the 
convective model always leads to positive Q\. 

In all our models, the convective parameters vary 
around the values defined in Table [2] As noted in Sec- 
tion l2.3l with such values we can reproduce the observational 
properties of the non-Blazhko RR Lyrae variables and also 
classical Cepheids (with smaller eddy- viscous dissipation). 
With these parameters we neglect the effects of turbulent 
pressure and overshooting from the convective zone. The 
computation of models including these effects is very time- 
consuming. In order to check whether the inclusion of tur- 
bulent pressure and overshooting can change our results we 
have computed one additional model sequence in which we 
set q p = 1.0 and at = 0.01. We adopted the same mod- 
ulation parameters as for set M6 (Tabled i.e., sinusoidal 
modulation with a period equal to 60 d and an amplitude 
of A = 50 per cent. The results are qualitatively the same 
as described in the previous sections. All models are charac- 
terised by positive values of Q\. For the lower temperatures, 
the light curve variation is qualitatively the same as pre- 
sented in, e.g., Fig. (set M6). At higher temperatures, the 
range of variation of the Fourier parameters becomes smaller 
for models including turbulent pressure. Consequently, for 
such models, an even larger amplitude of turbulent con- 
vection modulation would be necessary to reproduce the 
strongly modulated Blazhko variables. 



5 SUMMARY AND CONCLUSIONS 

In this pape r we have invest igated whether the mechanism 
proposed by IStothersl (^OOd 1 ) is capable of reproducing the 
light curve variation and the properties of the frequency 
spectra we observe in RR Lyrae Blazhko variables. The 
mechanism proposed by Stothers is extremely complicated, 
as it assumes the coupling between a variable magnetic field 
and turbulent convection in high amplitude, strongly non- 
linear pulsators. To get a variable magnetic field a rota- 
tional/turbulent dynamo is postulated. There are no com- 
prehensive models dealing with all these processes. Even the 
current models to deal with turbulent convection in non- 
linear pulsation are rather simple, ID, one-equation formu- 
las. 



Recent observational progress poses serious problems 
for the two most popular models to explain the Blazhko 
effect, the Magnetic Oblique Rotator/Pulsator model and 
Non-radial Resonant Rotator/Pulsator model (see Intro- 
duction). The Stothers mechanism remains as a scenario 
that has not been confronted with concrete challenges from 
observati ons. The variab le, tangled magnetic fields postu- 
lated by IStothersl ((2006) cannot be easily detected, mak- 
ing the idea hard to verify o n purely observational grounds 
l|Kolenberg fc Bagnuldbo'ogh . Its stochastic nature makes it 
attractive in the light of irregularities commonly detected in 
the Blazhko cycles of many variables. However, this is only 
a general idea, not supported by any detailed calculations 
or modelling. To advance with the theory we proposed a 
simple model to check whether the variation of convective 
structure of the star can lead to the modulation we observe 
in Blazhko stars. The modulation of turbulent convection is 
introduced into our models in an ad hoc way, neglecting the 
dynamical coupling with the postulated magnetic field. 

Comparison of our models with overall properties of 
Blazhko variables, as well as a detailed comparison with the 
strongly modulated prototype RR Lyr, observed by Kepler, 
reveals several discrepancies with the observations and chal- 
lenges for the Stothers mechanism. In our opinion, the most 
important objection is the required strength of the turbu- 
lent convection modulation. In order to reproduce the ranges 
of the light curve variation observed in strongly modulated 
stars like RR Lyr, we have to modulate the mixing-length 
by up to ±50 per cent on a time-scale of several tens of 
days. The same estimate results from th e analysis of pulsa- 
tion p eriod changes (Section l3.3l see also lMolnar fc Kollathl 
l2010h . The physical reality of such strong modulation is, in 
our opinion, questionable, although definite claims require 
detailed magnetohydrodynamic modelling of the problem, 
which is beyond the scope of this paper. We note here that 
it is not possible to infer the mixing-length variation in the 
subat mospheric l ayers directly from observations. In a recent 
paper, IPrestonl <|201lf) analysed the variation of FWHMs 
(full width at half maximum) of spectral lines in several 
RR Lyrae-type stars. This parameter is a measure of the at- 
mospheric turbulence. Preston shows that maximum value 
of FWHM varies with the Blazhko cycle. However, it is not 
clear if changing intensity of the atmospheric turbulence re- 
lates in any way to changes of the mixing-length in the sub- 
atmospheric layers. A strong variation of the FWHM oc- 
curs also in non-modulated RR Lyrae stars (see Preston's 
fig. 3), and, we may add, also in non- modulated hydro- 
dynamical models with c onstant mixing-length (see, e.g., 
iBenz fc Stellingwedll985h . In our opinion, the variation of 
the maximum FWHM reflects the behaviour of the turbu- 
lence generated by the velocity gradients in the atmosphere, 
rather than the possible variatio n of the mixing length in 
the deeper layers. IPrestonl (|201lT ) shows that the maximum 
value of the FWHM is strongly correlated with the pulsa- 
tion amplitude, which varies during the Blazhko cycle (his 
fig. 6). The larger the pulsation amplitude, the larger the 
maximum value of the FWHM. On the other hand, we note 
that the maximum of the FWHM occurs at pulsation phases 
close to the minimum radius. These are the phases at which 
the velocity gradients in the atmosphere are the strongest. 
Strong velocity gradients generate strong turbulence. Thus, 
the larger the pulsation amplitude, the stronger the atmo- 
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spheric velocity gradients and, consequently, the stronger 
the turbulence and the higher the FWHM maximumQ 

The light variation we compute resembles that observed 
in Blazhko variables, however, the details cannot be repro- 
duced, at least for RR Lyr. We note that our results can 
be used in the future for comparison with other Blazhko 
variables for which satellite data will be soon available. 
This would clarify how severe the discrepancies between the 
model and the observed light variation are. As for the fre- 
quency spectra, we cannot reproduce the asymmetry of the 
modulation side peaks around the main pulsation frequency. 
In our models, the higher frequency side peak always has a 
higher amplitude, which is not the case in a bout a quarter 
of the Blazhko variables (|Alcock et alj|2003h . 

The critical analysis of the Stothers idea by iKovacs! 

(2009) is al so worth mentioning. He points out that the pre- 
dictions bv lStoth ers (20061. 120101 ) of expected period changes 
through the instability strip, that agree with values observed 
in some Blazhko stars, are not correct, as they result from a 
misinterpretation of the linear and non-linear mode l periods. 

In the light of our results, the idea proposed bv lStothersl 

(2006) faces difficulties, and should be treated with caution. 
It needs refinement and further detailed studies in order to 
put it on solid physical grounds. Certainly, it is worth further 
investigation, but other possibilities to explain the Blazhko 
phenomenon should be explored as well. 
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